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Abstract. A method for exploring unstable structures generated by nonlinear 
dynamical systems is introduced. It is based on the sampling of initial conditions 
and parameters by Replica Exchange Monte Carlo (REM), and efficient both for the 
search of rare initial conditions and for the combined search of rare initial conditions 
and parameters. Examples discussed here include the sampling of unstable periodic 
orbits in chaos and search for the stable manifold of unstable fixed points, as well as 
construction of the global bifurcation diagram of a map. 



1. Introduction 

Nonlinear dynamical systems exhibit very complex structures in the phase space [H [2] • 
There can be a number of stable or unstable fixed points, periodic orbits, chaotic 
saddles, and more intricate invariant manifolds. Basin structures corresponding to 
these invariant manifolds can also be very complicated. Quest for these structures 
is an important subject in the study of dynamical systems. 

It is, however, not easy to capture these global structures by time evolution from 
randomly chosen initial conditions, because they are often unstable and correspond to 
very rare initial conditions. Analysis of these structures requires an efficient algorithm 
for sampling rare trajectories from rare initial conditions, which acts as an anatomical 
tool of nonlinear dynamical systems. Moreover, properties of chaotic or multi-basin 
dynamical systems are often sensitive to the choice of values of the parameters. Then, 
additional complexity arises when we do not have precise knowledge on the values of 
parameters with which interesting structures appear. Algorithm with which we can 
perform combined search of the space of initial conditions and of the parameter space 
is required for analyzing dynamical systems. 
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In this paper, we introduce a method based on Markov Chain Monte Carlo 
(MCMC) [21 m [5] , which is efficient both for the search of rare initial conditions and 
for the combined search of rare initial conditions and parameters. In the proposed 
algorithm, we identify initial conditions of trajectories as microscopic states to be 
sampled, and assign each of them an artificial "energy" that represents a measure 
of "atypicality" of the trajectory starting from the initial condition. These initial 
conditions are sampled by MCMC and a set of atypical trajectories is obtained. This 
approach is easily extended to the combined search of the space of initial conditions 
and of the parameter space. An essential idea is the use of direct product of initial 
conditions and parameter as a state space explored by MCMC. 

In both cases, however, naive applications of MCMC will be impaired by sensitive 
dependence on initial conditions and parameters in nonlinear dynamical systems, which 
lead to a highly multimodal energy function. Here we use Replica Exchange Monte Carlo 
(REM), which is also known as Parallel Tempering (PT), to circumvent the difficulty. 
REM is useful for the sampling on a rugged energy surface and intensively used for 
the finite-temperature simulation of spin glass [6], [7j and biomolecules [8|,[7]. When we 
generate multiple samples with conventional optimization methods, such as simulated 
annealing, it is difficult to control the probability of repeated appearance of the same 
(or similar) items in an obtained set of samples. The advantage of REM is that it 
realizes unbiased sampling from the given ensemble even in the region corresponding to 
low "temperatures". 

The idea of using MCMC or related algorithm for sampling rare trajectories in 
nonlinear dynamical systems is already seen in literature P [TOl [TTl [121 [131 [IH IIS] > while 
combined search of initial conditions and parameters is rarely studied. 

These studies are classified into two categories. Some of them [9l [TOl [Til [121 [13] . 
including an earlier work pj and recent studies [131 [12] , developed frameworks where the 
states to be sampled by the algorithm are entire trajectories or orbits, instead of initial 
conditions of the trajectories. The "transition path sampling" [161 IlHl [HI [H], which 
mostly used for sampling trajectories between two metastable states, is also based on 
a similar choice of state space, that is, a state consists of the array of positions (and 
momentums) of particles at all time-steps in a trajectory. The algorithm proposed in this 
paper is much simplified with the choice of initial condition as state variable. It is also a 
natural choice for a deterministic but not necessary reversible dynamical system. In the 
following sections, we will show that the proposed algorithm can deal with impressively 
various kinds of problems. 

Another type of algorithm proposed in [TH [15] approximates "pseudo trajectories" 
with desired property by a set of particles each of which obeys original equation of 
motion. Using a genetic algorithm like split/delete procedure, efficient sampling of 
rare trajectories is realized. It is a kind of sequential Monte Carlo [HI [19] or diffusion 
Monte Carlo [201 HI] and not genuine MCMC. It can also be interpreted as a multi- 
particle version of ingenious "staggered-step" algorithm [21] developed earlier. The 
idea seems effective for trajectories with positive Lyapunov numbers, but its advantage 
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might be reduced in a search for trajectories with negative Lyapunov numbers where 
spht trajectories does not diverge fast without strong external noise. Anyway, it seems 
to capture different aspects of chaotic systems and will be complementary to the method 
proposed in this paper. 

The organization of the paper is as follows: In Section. 2, we explain the basics 
of the proposed method. In Section. 3, two examples of initial condition sampling are 
discussed. The first one is a toy example of the sampling of the unstable periodic 
orbits of the Lorenz model. The second one is sampling of stable manifold of the 
unstable fixed points of a double pendulum with dissipation, which is a much difficult 
and interesting example. In Section. 4, search of parameter space and combined search 
of parameter space and initial condition space with a modified algorithm are studied. 
Section. 5 devoted to the discussion on the results and possibility of further extensions. 

2. Method 

2.1. State Vector 

In the basic algorithm for the discrete time dynamics with a given set of parameters, 
the initial condition X = {xi,...,Xn) G M" of a trajectory is used as a state variable 
that characterizes the trajectory evolving from the initial condition. In a continuous 
time case, it is often better to include the time T when desired phenomena occur. Then 
a state will be X = {xi, . . . ,Xn,T) G M"^^. Finally, when the parameters are also 
unknown, the state vector will be X = (xi, . . . , a;„, T, ai, . . . , a^) G K^+i+^-j where 
a = («!, . . . , am) is a vector of unknown parameters. 

2.2. Energy Function and Gibbs Distribution 

In order to explore rare structures in phase space, we should construct a fictitious 
"energy" function E{X) defined on the state space. The function depends on the kind of 
rare orbit we want to explore. For the detection of periodic orbits of a map Xn+i = f{xn), 
a simply possibility of the energy function is E{xo) = log(|/"'(a:o) — /(xo)|), where the 
state X coincides with the initial condition xq of the iteration. The other choices of the 
energy function to explore atypical structure in phase space are shown in Sec. HI where 
E may depend also on T and a. 

Once we define the energy function, it is straightforward to define the Gibbs 
distribution with the energy 



While p{X) coincides with the uniform density in a prescribed region when /3 takes a 
small value, it concentrate to regions with small values of E when jS becomes large, 
where we can collect samples of trajectories (and parameters) of desired properties. 



p{X) 



exp{~(3E{X)) 

m 




(1) 
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2.3. Metropolis Algorithm 

Now, the problem of sampling atypical trajectories is reduced to the sampling from the 
density p{X). Here we use the Metropolis algorithm [22], a simplest implementation of 
the idea of MCMC The Metropolis algorithm used here is essentially the same as the 
standard one commonly used in statistical physics. There is, however, an important 
difference in actual implementation, that is, we should simulate the trajectories of the 
length T from the initial condition x with the parameter a at each trial of changing x, T, 
and a. If we explicitly represent this difference, the iteration of algorithm is described 
as follows. 

(i) Draw a perturbation to current states (Ax, AT, Aa) from a prescribed "trial" 
density q, which defines a move set. Hereafter, the mirror symmetry 

q{-Ax, - AT, - Aa) = q{Ax, AT, Aa) 

of the density q is assumed. Then, when the current values of (x, T, a) is 
(x*-"-*, T*^"\ a*-"-*), the candidate of the next states is defined as {x',T',a') = 
(x(") + Ax, T(") + AT, a(") + Aa). 

(ii) Run the simulation of length T' of the dynamical system with parameter a' from 
the initial condition x'. From the result, E{a',T',x') is computed. 

(iii) Draw a uniform random number R G [0, 1]. If and only if 

exp(-/3E(x^T^aO) 
exp(-/5E(xW,TW,aW)) ' 

the new proposal is accepted: x*^""''^'' = x', T*^""*"^^ = T', a*^""^"*^^ = a'. Else nothing is 
changed: x('"+^) = x^"'\ T("+^) = T^"), a("+^) = a^") 

(iv) Return to Step (i). 

An important point for treating unstable structures in potentially chaotic systems is 
the choice of the density q, or, equivalently, the set of moves. The point is that the moves 
should be hierarchical for the efficient sampling, that is, the coexistence of tiny and large 
changes in phase space is essential [21]. Here we adopt the method introduced in [21] . 
that is, the elements Axj of the perturbation Ax is given by 6xi = s x dx 10~^ where d 
and e are random integers uniformly distributed in d e [1, 9] and e G [A^™°, N^^^], and 
s is a binary random number that takes the value ±1 with probability 0.5, respectively. 
The parameters A*"™™ and N^^^ of the algorithm define the logarithmic scales of the 
largest and smallest perturbations. The corresponding trial density g(xj) becomes a 
mixture of uniform distributions with the different scales. It has a sharp peak near zero 
as well as a very long tail. The components of Aa and AT are also generated by a 
similar way. 

This version of Metropolis algorithm appears to provide a simple and universal way 
of treating the Gibbs distribution ([1]) . The efficiency of the algorithm, however, can be 
reduced when /3 become large in the case of a highly multimodal energy function. This 
difficulty will be treated by replica exchange Monte Carlo algorithm described in the 
next subsection. 
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2.4- Replica Exchange Monte Carlo 

Replica exchange Monte Carlo (REM) provides an efficient way to investigate systems 
with rugged free-energy landscapes [23l EH E] , particularly at low temperatures. In the 
present context, it is used in references [12], [13] for the sampling of unstable orbits. It 
is also used in the framework of "transition path sampling" in references [TTl [T7] . 

In a replica exchange Monte Carlo simulation, a number of systems {X(m)} with 
different inverse temperatures (3rn {replicas) are simulated in parallel. At regular 
intervals, an attempt is made to exchange the configurations of selected, usually 
adjacent, pair of replicas. It is accepted with the probability 

P(X(„+i) ^ X(^)) = min[l,exp(A/5AE)], 

where A/5 = Pm+i — Pm is the difference between the inverse temperature of the two 
replicas and AE = E{X(^rn+i)) — -£'(-^(m)) is the energy difference of them. 

The exchange of replicas with different temperatures effectively reproduces repeated 
heating and annealing, which avoids trapping in local minima of the energy. Note that 
it is especially useful with hierarchical moves defined in the previous subsection, because 
large moves are accepted at high temperatures and tiny moves are dominated at low 
temperatures. 

On the other hand, the above rule of stochastic exchange preserves the joint 
probability distribution 

m 

as shown in the literature, e.g., ^J. Thus, even with the annealing effect, the 
probability distribution of each replica coincides with P^,^(X(m)), which means 

that an unbiased set of samples is obtained at all inverse temperatures {Pm}- 

3. Initial Condition Sampling 

In this section, we give examples of the initial state sampling by the proposed method. 
An example is the search for unstable periodic orbits (UPOs) of the Lorenz model. 
Another interesting example is sampling of rare orbits in a double pendulum system, 
i.e., initial conditions that locate on the stable manifold of unstable fixed points are 
sampled by the proposed method. 

3.1. Unstable Periodic Orbits of the Lorenz model 

In this subsection, we show that unstable periodic orbits (UPOs) in a continuous-time 
dynamical system can be detected by the proposed method. UPOs are considered 
as important objects that govern the properties of chaos [2] and there are considerably 
many works that deal with the computation of UPO [25], EH] ET] EH [121 IB] • Our purpose 
here is to show that how the proposed method works with a familiar example, but not 
to prove that the proposed method is superior to all of these existing methods. It will 
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be, however, useful to note that the proposed method with REM can generate UPOs 
uniformly under the assumption of uniform measure in the space of initial conditions. 
This suggests that it can be useful for the global search for UPOs in combination with 
other methods. 

Here, we consider the Lorenz equations ^291 [30] 

x = a{y- z) 

y = rx — y — xz (2) 
z = xy — bz, 

as a typical continuous-time chaotic system, where a, r and b are parameters of the 
system. The state of the system is denoted in a vector format as x = (x, y, z) and the 
initial condition is written as xq = {xo,yo,zo). The flow generated by the dynamical 
system ([2]) is expressed as and the orbit determined by the initial condition xq is 
written as Xt = (p^{xo). 

The state sampled by REM is defined by X = {xo,T) = {xo,yo, zo,T). Then, a 
candidate of the energy function is given as 

E{xo,T) = \og{\^^{xo)-xo\+e). (3) 

The parameter e <^ 1 is a constant for avoiding the divergence of the energy. When the 
orbit {ip^{xq) : t G M} is periodic, there exists a time T such that (p'^{xo) = xo, and the 
energy of the initial condition xq takes the minimum value log(e). 

There are, however, problems with the naive choice ([3]) of the energy function. First, 
the energy E{xo, 0) always takes the minimum value, because (f^{x) = x. Also, when an 
initial state Xq locates on a fixed point, the energy of the initial state takes the minimum 
value for all T. This implies that the almost all initial conditions sampled by the above 
energy function will be in the vicinity of fixed points. To avoid these unfavorable 
situations, we will add a penalty term P{xo,T) to the naive energy function ([3]). An 
improved energy function is 

E{xo,T)= log(|<^^(xo)-fo|+i^(^o,T) + e) (4) 
P{x,,T)=g(^ V-{xofdT,v,^+g{T,T,), (5) 

where we use an auxihary function g{s,Sc) = 6(sc — s){^ — j-). 9(s) is the Heviside 
function defined by 0(s) = 1 if s > else 0(s) = 0. The first term in the equation ([5]) 
represents a penalty for slower "average speed" of the trajectory, where Vc is a threshold 
parameter. When an initial condition is in the vicinity of a stable or unstable fixed 
point, the orbit stays near the point within a certain time. Then, the averaged speed 
becomes slower and the value of integral in the penalty term becomes smaller, which 
causes a large value of the energy. The second term is a penalty to the states that have 
small T, where Tc is a threshold parameter. 

For sampling of initial condition xq of the Lorenz model, we take account of the 
symmetry {x,y) {—x,—y) of the equations. We sample initial conditions from the 
Poincare section xq = {xo,yo,r — 1), where Xq G [0,20] and yo G [0,20]. The period 
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Figure 1. (a) Samples generated by the proposed method are plotted on the {E,T) 
plane. The crosses, triangles, filled circles, and circles correspond to /? = 2,3,4,5, 
respectively, (b) Samples with the energy lower than —1.5 are plotted on the Poincare 
section (x,y,r — 1). The clusters, each of which corresponds to different minima, are 
identified by a cluster analysis method and they are plotted by a different symbols. 25 
clusters are shown. The inset is a magnification, e = 10^^, A^™™ — 1, and iV™"^ = 5. 

T of orbits is assumed to be in the interval [0,6]. Using the energy function ([5]), 
the states X = {xo,T) are sampled by the proposed method. 31 replicas with 
(3 = 2.0 + O.li, i = 0, 1, . . . , 30 are used. 

The samples obtained by the proposed method are plotted on the {E, T) plane in 
Fig. [l](a) for different temperatures, which are obtained in parallel in a replica exchange 
Monte Carlo calculation. For lower temperatures, the energy takes smaller values for 
specific values of the period T, which correspond to UPOs. 

In Figure[l](b), the samples of initial conditions with E{xo, T) < —1.5 are plotted on 
the Poincare section. The initial conditions form clusters and each cluster corresponds 
to initial condition in the vicinity of a different UPO. In each cluster, we pick up the 
initial condition that has the minimum energy and calculate orbit ip^{xo) starting from 
the initial condition. Typical orbits are shown in Fig. O These orbits are closed in high 
precision, i.e., the difference between xq and (p'^{xo) is in order of 10~^, indicating that 
they are UPOs. 

By using AUTO software[3l], we tested these orbits. An orbit detected by the 
proposed algorithm is used as an initial guess of a Newton's method implemented in 
AUTO. We see that the convergence is remarkably quick, which indicates that the 
iteration begins "sufficiently near" UPOs. The difference between the output of the 
Newton's method and the initial guess is also very small. 




Figure 2. Typical unstable periodic orbits of the Lorenz model sampled by REM are 
shown. Each orbit is produced from time evolution starting from the initial condition 
that has the minimum energy in a cluster shown in Fig. [IJb) . 



3.2. Stable Manifold of Unstable Fixed Points of a Double Pendulum 

Let us consider the following dynamics of the double pendulum with dissipation 

_ sin (2 (gi - e^)) + 2 sin {9^ - O^) 9^ + 3 sin (gQ + sin {Oi - 262) _ ■ 

cos (2 (01-^^2)) -3 

2 sin (^1 - 62) (261 + cos (^1 - 62) 62 + 2 cos I 



= cos (2 (g,-^,))- 3 

where A; is a dumping coefficient. 

The double pendulum system ([6]) has three unstable fixed points: {61,61,62,02) = 
(7r,0,7r,0), (tt, 0,0,0) and (0,0,7r,0). Each unstable fixed point corresponds to an 
"inverted" state of pendulums. Starting from any initial condition xq = {61,61,62,62) = 
{0,uJi,0,uj2), however, almost all trajectories are converged to the stable fixed point 
{61, 61, 62, 62) = (0, 0, 0, 0) by dissipation originate from the friction at the hinge. 

In this system, we try to detect atypical trajectories that converge to unstable 
fixed points. When an initial condition locates on the stable manifold of an unstable 
fixed point, the orbit starting from such an initial condition converges to the unstable 
fixed point after a long time evolution. We search such an atypical initial state on the 
Poincare section {61,61,62,62) = (0, tui, 0, u;2)- 
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The state sampled by the proposed method is set as X = {u!i,uj2,T). There are 
many possibihty of the energy function. Here we try to find trajectories converging one 
of the three fixed points and choose the energy as 



where Ej, represents "kinetic energy" of the pendulum at time T stating from an initial 
condition defined by (ci;i,ci;2)- The term Ep represents an artificial "potential energy", 
which has minimum at each of the three unstable fixed points. Note that if we want to 
find trajectories converging to a given unstable fixed point, we can use other artificial 
potential energies such as i?p(to'o, cui, T) = |6'i(T) — 7r| + |6'2(r) — 7r|. The energy landscape 
defined by ([7]) is shown in Fig. [3](a), where the energy E{u!i,u!2,T) vs. ui and uj2 is 
plotted for a fixed time, T = 5. Rugged structure is clearly seen and increase of the 
sampling efficiency by REM is expected. 

We discuss the result of a simulation with 31 replicas with /3j = 2.0 + 0.2i for 
i = 0, 1, . . . , 30. In Fig. Mjo), initial states sampled by the proposed method is plotted 
on the {uji, E) plane for replicas with inverse temperatures (3 = 0.1, 1.1, 2.1, 3.1, 4.1, 5.1. 
The points are scattered at higher temperatures, but they are concentrated in separated 
regions when the temperature becomes lower. 

Because states with lower energies correspond to desired atypical trajectories, we 
select a set of the states whose energies E are lower than a threshold E^r = —5. We show 
these states in Fig. Hl^a). By using cluster analysis, we divide these initial states into 
clusters. Since these clusters are well separated in the plane, we identify each cluster as 
a representation of a qualitatively different set of trajectories. 

In each cluster, we identify an initial state that has minimum energy. Figure 111(b) 
demonstrates atypical trajectories starting from these initial conditions on (^i, 62) plane. 
It is seen that each cluster corresponds to qualitatively different pattern of motion. 
Examples of sequences of snapshots for pattern of the double pendulum are shown in 
Fig. 111(c), each of which corresponds to a trajectory obtained by the above procedure. 

4. Parameter Search and Combined Search 

In this section, search in a parameter space and combined search in parameter and initial 
condition spaces are discussed as an extension of our method. Examples are sampling of 
the boundary of the Mandelbrot set and exploration of the global bifurcation structure 
of the logistic map. 

4.1. The Boundary of the Mandelbrot Set 



For a map fc{z) = + c parameterized by a parameter c G C, consider the sequence 
(0, /c(0), fdfci^)), • • •)■ The Mandelbrot set [32] is defined as the set of points c G C 



E{uJi,uj2,T) 

Ek{uJi,uj2,T) 

Ep{uJi,uj2,T) 



log{Ek{LJi, UJ2, T) + Ep{uJi, ^2, T)), 

\e,{T)\ + \el{T)i 

min[ cos(^i(T)) cos(^2(T)), cos(^i(T)) + cos(e2(T)) ] + 1, 



(7) 
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Figure 3. (a) Rugged energy landscape E{xq,T) with fixed T — 5.0 as a function 
of uJi and lu2. (b) Samples are plotted on {E,uji) for different temperatures; (3 = 
0.1, i.l, 2.1, 3.1, 4.1, 5.1. The sampling is performed by 51 replicas with (3 = O.li for 
i = 1, 2, . . . , 51. 10-^ TV™" = 1, and TV™*^^ = 9. 



such that the above sequence does not escape to the infinity. Hereafter, we denote 
WM) = fm) and /c(/c(/c(0))) = mO), and so on. 

The Mandelbrot set is seen to have an elaborate boundary in the complex plane, 
which does not simplify at any magnification. This qualifies the boundary as a fractal. 
In order to compute the boundary of the set, we use the parameter c as a state of the 
proposed method and employ the following energy function 

E{c) = -\og{n{c)), 

where the function n{c) is defined as the smallest n that |/"(0)| > 2. If n > N for a 
prescribed large number = 200, we set n(c) = 1 and E{c) = 0. In that case, the point 
c is considered as inside the Mandelbrot set and not on the boundary. 

Using this energy function, the boundary of the Mandelbrot set is calculated by 
the proposed method and the results are shown in Fig. [5l While the points are almost 
randomly distributed in the complex plane for the replicas with a high temperature 
l/jS, the distribution of points corresponding to replicas with lower temperature is 
concentrated on the boundary of Mandelbrot set. 

4-2. Periodic Orbits and Bifurcation Diagram of the Logistic Map 

The logistic map [331 [I] is defined by 

fa{x) = ax{l - x), (8) 
where a G [0, 4] is a parameter. Let us consider the energy function 

E{a,xo) = \og{\f^{xo) -xo\+e), (9) 

where A; is a given period and e is a constant that determines the minimum energy. 
Indeed, we can find initial conditions xq corresponding to period k orbits for a given 
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Figure 4. (a) States with i? < — 5 is plotted on tlie (wi,W2) plane. The states form 
clusters, which are regarded as sets of qualitatively different trajectories. The different 
cluster is plotted with different symbols in the figure, (b) Typical trajectories starting 
from the state that has minimum energy within a cluster are plotted on the {9 1,62) 
plane. Six trajectories are shown, (c) Sequences of snapshots of pattern of the double 
pendulum that correspond to trajectories. Three trajectories, each of which converges 
to a different fixed point are shown in the upper, middle, and lower panel. The other 
parameters are the same as those in Fig. [3] 



parameter a using the energy function ([9]) by the proposed method (the resuhs are 
not shown here). On the other hand, the extension of samphng state space from xq 
to (a, xq) enables us the study of global bifurcation structure of the periodic solutions. 
The extension is straightforward, i.e., we sample the vector (a, xq) from the Gibbs 
distribution determined by the same energy function ([9]) using REM. 

In Fig. [6l the result for the period = 3 is shown, where sampled points are plotted 
on the (a, xo) plane. We sample orbits from a G [3.8,4] and x G [0, 1]. The points are 
scattered broadly in the vicinity of the true bifurcation branch of the period /c = 3 
orbits at higher temperatures. The dispersion of the points becomes smaller at lower 
temperatures and periodic orbits with = 3 are detected. 

5. Summary and Discussion 

In this paper, we presented a general strategy for exploring unstable structures generated 
by nonlinear dynamical systems. An artificial "energy" is defined as a function of the 
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(b) 



0=1- 




Re(c) 



Figure 5. (a) The Mandelbrot set. The black region shows the Mandelbrot set and 
the outside of the set is colored according to the escape time, (b) The boundary of 
the Mandelbrot set. Six replicas with (3 — 1.0,2.0,3.0,4.0,5.0, and 6.0 are used for 
sampling. TV^""' = 1 and TV™*^^ = 9. 




Figure 6. The bifurcation structure of the period three solutions of the logistic 
map is obtained by the proposed method. Ten replicas are used with P = 0.5i for 
i = 1, 2, . . . , 10. e = 10-*, TV™" = 1, and iV™"^ = 11. (a) (3 = 0.5. (b) (3 = 1.0. 

initial condition of the trajectory and the Gibbs distribution induced by the energy 
is sampled by the Metropolis method. When an "energy" suitable for the purpose 
is chosen, sampling from the Gibbs distribution at low temperatures realizes efficient 
sampling of rare initial conditions that leads to interesting trajectories. Replica exchange 
Monte Carlo (REM) is used so as to avoid capturing in local minima of the energy 
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landscape. 

While implementation of the proposed method is simpler than the methods in which 
artificial energy is defined in the space of entire trajectories, the method can treat a 
variety of problem, including the search for unstable periodic orbits of the Lorenz model 
and stable manifold of unstable fixed points of a double pendulum. 

It is also shown that search in the parameter space and combined search of initial 
conditions and parameters become possible by adding the parameters to the state vector. 
Two examples shown here are the Mandelbrot set of a complex dynamical system and 
the bifurcation diagram of the logistic map. 

An important future problem is to calculate quantitative results by the proposed 
method. Examples are the relative densities of initial conditions that lead to different 
Lyapunov numbers and densities of escape time from a chaotic saddle. These 
calculations are possible because REM is not only an optimization method, but also 
realize correct sampling from the Gibbs distribution. It is also interesting to introduce 
other types of extended ensemble Monte Carlo, such as multicanonical algorithm, for this 
purpose. Research in this direction is now in progress and the results will be published 
elsewhere. 

The combined search of initial conditions and parameters by the proposed method 
is also promising and should be tested in the study of more complicated systems. There 
is, however, an inherent interpretation problem, i.e., the joint density in the direct 
product of parameter space and initial condition space seems not to have a definite 
interpretation. The method is already useful to give a rough sketch of the bifurcation 
diagram, but it will be better if we can give a physical meaning to the joint density. 
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